#ifndef __STVENANTKIRCHHOFF_H
#define __STVENANTKIRCHHOFF_H 1

#include <feel/options.hpp>
#include <feel/feelcore/application.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.hpp>
#include <feel/feelpoly/im.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelmesh/meshmover.hpp>
#include <feel/feelvf/vf.hpp>
#include <feel/feeldiscr/bdf2.hpp>
#include <feel/feeldiscr/operatorinterpolation.hpp>
#include <research/life-workspace/Fluid-Structure/functionSup.cpp>

namespace Feel
{

template<int Dim,
         int Order,
         template<uint16_type,uint16_type,uint16_type> class Entity>
class StVenantKirchhoff : public boost::enable_shared_from_this< StVenantKirchhoff<Dim, Order, Entity> >
{
public:

    typedef StVenantKirchhoff<Dim, Order, Entity> self_type;
    typedef boost::shared_ptr<self_type> self_ptrtype;

    typedef double value_type;

    typedef Application application_type;
    typedef boost::shared_ptr<application_type> application_ptrtype;

    typedef Backend<value_type> backend_type;
    typedef boost::shared_ptr<backend_type> backend_ptrtype;

    /*matrix*/
    typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    typedef typename backend_type::vector_ptrtype vector_ptrtype;

    /*mesh*/
    typedef Entity<Dim,1,Dim> convex_type;
    typedef Mesh<convex_type> mesh_type;
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

    /*basis*/
    typedef Lagrange<Order, Vectorial> basis_u_type;
    typedef bases<basis_u_type> basis_type;

    /*space*/
    typedef FunctionSpace<mesh_type, basis_type,Continuous, double> space_struct_type;
    typedef boost::shared_ptr<space_struct_type> space_struct_ptrtype;

    BOOST_MPL_ASSERT( ( boost::is_same<typename space_struct_type::bases_list, basis_type> ) );

    /* functions */
    typedef typename space_struct_type::element_type element_struct_type;
    typedef boost::shared_ptr<element_struct_type> element_struct_ptrtype;

    typedef element_struct_type element_dirichlet_type;

    typedef Bdf<space_struct_type>  bdf_struct_type;
    typedef boost::shared_ptr<bdf_struct_type> bdf_struct_ptrtype;

    /* export */
    typedef Exporter<mesh_type> exporter_type;
    typedef boost::shared_ptr<exporter_type> exporter_ptrtype;

    StVenantKirchhoff( application_ptrtype __appli );

    application_ptrtype application() { return M_lifeApplication; }

    void changeRepository(std::string pathLoc="");

    mesh_ptrtype createMesh( double meshSize );

    //-----------------------------------------------------------------------------------//

    void init();

    void solveHyperElasticity();

    void updateCLDirichlet(vector_ptrtype& U);

    //void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
    void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J, vector_ptrtype& R );
    void updateResidual( const vector_ptrtype& X, vector_ptrtype& R );

    void updateTime(double __t);


    void updateVelocity();

    template <typename element_fluid_ptrtype>
    void updateStressTensor(element_fluid_ptrtype fluidSol);

    template <typename element_fluid_ptrtype>
    void updateVelocityInterface(element_fluid_ptrtype fluidSol);


    void updateBdf();

    void exportResults();

    void solveElasticity();
    void updateLinearPDE(const vector_ptrtype& X,sparse_matrix_ptrtype& A , vector_ptrtype& F);


    void predictorDispl();

    //-----------------------------------------------------------------------------------//

    boost::shared_ptr<element_dirichlet_type>  getDirichletElement();

    space_struct_ptrtype functionSpace() { return M_Xh; }

    backend_ptrtype backend() { return M_backend; }

    element_struct_ptrtype getSolution() { return U_displ_struct; }

    element_struct_ptrtype getVelocity() { return U_velocity_struct; }



    element_struct_type & meshVelocity() { return *(getDirichletElement()); }

    bool weakCL() { return M_weakCL; }

private:


    //void solve( sparse_matrix_ptrtype const& D, element_struct_type& u, vector_ptrtype const& F, bool is_sym );

    //void exportResults( element_struct_type& u, element_struct_type& v );
    void exportResultsInstat( element_struct_type& u, element_struct_type& v,double time );

private:

    application_ptrtype M_lifeApplication;

    backend_ptrtype M_backend;

    double meshSize;
    mesh_ptrtype M_mesh_struct;

    MeshMover<mesh_type> M_mesh_mover;

    space_struct_ptrtype M_Xh;

    element_struct_ptrtype U_displ_struct;
    element_struct_ptrtype U_velocity_struct;
    element_struct_ptrtype U_acceleration_struct;

    element_struct_ptrtype M_velocityInterface;

    bdf_struct_ptrtype M_bdf_displ_struct;
    bdf_struct_ptrtype M_bdf_velocity_struct;
    bdf_struct_ptrtype M_bdf_acceleration_struct ;


    element_struct_ptrtype M_normalStressBoundary;

    bool M_weakCL;
    //coef de penalisation
    double penalbc;

    double M_time;

    //Parametre du modele
    double EE;
    double sigma;
    double mu;
    double lambda;
    double M_rho;

    //Exporteur
    exporter_ptrtype M_exporter;
    uint M_cpt_export;

    bool M_LookComment;

}; // StVenantKirchhoff


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    template <typename element_fluid_ptrtype>
    void
    StVenantKirchhoff<Dim, Order, Entity>::updateVelocityInterface(element_fluid_ptrtype fluidSol)
    {
        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateVelocityInterface start\n";

        this->changeRepository();

        auto ola = fluidSol->template element<0>();

        auto dirElem = this->getDirichletElement();
        modifVec222(*dirElem,M_velocityInterface, vf::idv(ola) , "Interieur" );

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateVelocityInterface finish\n";

    }



    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    template <typename element_fluid_ptrtype>
    void
    StVenantKirchhoff<Dim, Order, Entity>::updateStressTensor(element_fluid_ptrtype fluidSol)
    {

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateStressTensor start\n";

        this->changeRepository();

        auto ola = fluidSol->template element<0>();

        auto dirElem = this->getDirichletElement();
        modifVec222(*dirElem,M_normalStressBoundary, vf::idv(ola) , "Interieur" );

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateStressTensor finish\n";

    }



} // Feel
#endif /* __NSSTOKES_H */
